Nonperturbative Approach to Circuit Quantum Electrodynamics 
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We outline a rigorous method which can be used to solve the many-body Schrodinger equation for 
a Coulomb interacting electronic system in an external classical magnetic field as well as a quantized 
electromagnetic field. Effects of the geometry of the electronic system as well as the polarization 
of the quantized electromagnetic field are explicitly taken into account. We accomplish this by 
performing repeated truncations of many-body spaces in order to keep the size of the many particle 
basis on a manageable level. The electron-electron and electron-photon interactions are treated in 
a nonperturbative manner using "exact numerical diagonalization". Our results demonstrate that 
including the diamagnetic term in the photon-electron interaction Hamiltonian drastically improves 
numerical convergence. Additionally, convergence with respect to the number of photon states in 
the joint photon-electron Fock space basis is fast. However, the convergence with respect to the 
number of electronic states is slow and is the main bottleneck in calculations. 
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I. INTRODUCTION 

To describe the interaction between matter and a 
single-mode quantized electromagnetic field, some ver- 
sion of the Jaynes-Cummings (JC) model is often ap- 
plied pp. The JC-model was first employed by Jaynes 
and Cummings to describe the interaction of photons 
with molecules but since then it has also been used in 
cavity electrodynamics to successfully describe matter- 
photon interaction in semiconductor nanostructures such 
as quantum dots [2J and in superconducting qubits [3l |4] . 
Advances in the field of circuit quantum electrodynamics 
have enabled us to enter the ultrastrong coupling regime 
where the photon-matter coupling strength reaches a con- 
siderable fraction of the energy of a single cavity pho- 
ton. This has been achieved by taking advantage of large 
dipole moments and long coherence times in supercon- 
ducting flux qubits (SHZ] and semiconductor quantum 
wells [MTO] embedded in high quality micro-cavities. 

In the ultrastrong regime, the JC model fails and evi- 
dence of the breakdown of the JC-model with the rotating 
wave approximation has been observed experimentally 
in superconducting [Hj and semiconductor systems [HI [P] . 
Exact numerical calculations predict the failure of the 
JC-model (even without the rotating wave approxima- 
tion) where the effects of the diamagnetic matter-photon 
interaction term as well as effects of states which are not 
part of the two level system approximation come into 
play with high coupling strength [11 . 

Using the method described later in this publication. 
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we have been able to calculate time dependent electron 
transport through a photon cavity [12] and to test the va- 
lidity of the Jaynes-Cummings model in the ultrastrong 
coupling regime |11| . With our approach, it would be 
relatively easy to add a time dependent perturbation to 
the closed system and integrate the equation of motion 
numerically. Choosing the frequency of the perturbation 
such that the EM field does not have time to adjust adi- 
abatically, it is possible to investigate non-adiabatic dy- 
namics related to the dynamical Casimir effect [T3] where 
photons can then be excited out of vacuum in correlated 
pairs. This non-adiabatic effect was recently observed 
experimentally for the first time [14]. 

In this paper we describe a general method which 
can be used to describe the interaction between an elec- 
tronic/atomic system with a single-mode quantized elec- 
tromagnetic field. We begin by calculating eigenfunc- 
tions and energies of the single-electron Hamiltonian (ini- 
tially completely ignoring many-body effects and the EM 
field) . We then use a number of the lowest single-electron 
eigenstates to construct a many-electron Fock state basis 
which is used to compute the eigenstates and energies of 
the many-electron Hamiltonian including the Coulomb 
interaction between electrons. Finally, we use a num- 
ber of the lowest Coulomb interacting eigenstates to con- 
struct a joint electron-photon basis. In diagonalizing the 
electron-photon Hamiltonian we obtain its eigenstates 
which include the electron-photon and electron-electron 
interaction "exactly" in the sense that the only approxi- 
mations are the finite sizes of single/many particle bases 
and finite size of grids used for numerical integration. 
The results are convergent with respect to these param- 
eters in a controllable manner. 

The paper is organized as follows. In Sec. ITT] we give a 
description of the single-electron Hamiltonian and cal- 



culate its eigenfunctions, which we use as a basis for 
many-body calculations. In Sec. 



Ill we introduce the 



second quantization many-body formalism needed to ac- 
count for the Coulomb interaction between electrons. In 
Sec. |IV| we couple the electronic system to single-mode 
quantized electromagnetic field and solve the many-body 
Schrodinger equation using a basis of Coulomb interact- 
ing electron states as well as photon Fock states. Results 
and concluding remarks are presented in Sees. [V|and [VT| 
respectively. 



II. SINGLE-ELECTRON HAMILTONIAN 

The system under investigation is a two-dimensional 
electronic nanostructure exposed to a static (classical) 
external magnetic field at a low temperature. The elec- 
tronic nanostructure is assumed to be fabricated by split- 
gate configuration in the y-direction, forming a parabolic 
confinement with the characteristic frequency f2o on top 
of a semiconductor heterostructure. The ends of the 
nanostructure in the x-direction at x — ±Lx/2 are 
etched, forming a hard- wall confinement of length Lj.. 
The external classical magnetic field is given by B = Bz 
with a vector potential A = {—By, 0, 0). Since we are in- 
terested in geometrical effects, we need the single-electron 
eigenstates to construct a many-body basis. We therefore 
need to solve the time independent Schrodinger equation 
for the Hamiltonian 
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where m is the effective mass of an electron 
charge, p the canonical momentum operator, ujc 
is the cyclotron frequency and ili^, = y'uj^ 
modified parabolic confinement. Note that the spin de 
gree of freedom is neglected. With the boundary con- 
ditions 'ip{±Lx/2,y) — ip{x,±oo) — 0, the mixing term 
iujcypx makes it impossible to use separation of variables 
to solve the time independent Schrodinger equation for 
the Hamiltonian in Eq. (Ill analytically. This means we 
will have to resort to numerical techniques. This proce- 
dure is relatively straightforward and will only be briefly 
covered here. 

To solve the time independent Schrodinger equation for 
-ffo, we compute the matrix representation of Hq in the 
basis {\(t>n) ® Wm)} where |(/()„) (g) \ipm) are eigenstates of 
Hq when the mixing term iuicyPx is omitted. The matrix 
elements are calculated analytically. Furthermore let us 
assume we have a bijection (n, m) — > i such that we can 
label the basis states using a single index i such that 
\^i) ~ \<t>ni) ® Wriii)- In coordinate representation, we 
have 



{x\(j)n^) = 




if rij = 1,3,5, ... 
if n, = 2,4,6,... 
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where a„ — ^Jhjm^l^ is a characteristic length of the 
system and _ff,„. are Hermite polynomials. 

After computing the matrix representation of Bq in 
the chosen basis, we diagonalize it and obtain it's eigen- 
states 1-0^) and corresponding energies Ei which satisfy 
HoVi^i) = ^ilV'i)- Note that 1-01) is the ground state, 102) 
the first excited state etc. In the diagonalization process 
we also obtain a unitary transformation which satisfies 
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Finally the wavefunctions of the lowest A'sos single- 
electron states '0i(r) are calculated and saved on a grid 
using 
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where iV^jj, is the number of basis states used for calcula- 
tions. In actual calculations we used approximately 120 
basis states in the x-direction and 31 in the y-direction 
so n e [1, 120] and to G [0,30], giving N^y = 120 x 31 = 
3720. This is a large enough basis such that numerical 
error due to the truncation is much smaller than the er- 
ror due to later truncation of many-body spaces. For 
this reason we will not investigate convergence for the 
single-electron system in this paper. 



III. MANY-ELECTRON HAMILTONIAN 

We can write the many-electron Hamiltonian as a sum 
of two terms 'He = H^ -I- "He where "He only contains 
the Coulomb interaction between electrons. Using the 
single-electron eigenstates |0i) = |i) as a basis, we can 
write the two terms in second quantization as |15| 

ni = Y,mo\j)4dj - J2 ^^4d^ (6) 

ij i 

He = \j2^t3\Vc\rs)dld]d,dr (7) 

ijrs 

where d} (di) are fermionic creation (annihilation) oper- 
ators of an electron in state \i) . The operators satisfy the 
usual fermionic anti-commutation relation {di,d-} = Sij 
and all other anti-commutators are zero. The matrix el- 
ement {ij\Vc\rs) in Q is a double integral in the spacial 
variables and involves integration with respect to the ob- 
servation location r 



{ij\V\rs) 



dr0*(r)J,,(r)^,(r) 
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and the integration with respect to the source location r' 



X,>(r)= dr'rAr')Vc{r,r')Mr') 



where Vc is the Coulomb potential given by 



Vc{r,r') 
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where 77 is a small positive regularization parameter. The 
integrals in (Is]) and ^ can not be done analytically due 
to the nontrivial geometry so they are performed numer- 
ically using a Gaussian quadrature scheme. We have to 
be careful with the numerical integration because techni- 
cally the wave functions reach infinity in the y-direction. 
although exponentially decaying. We therefore have to 
find some sensible cutoff in the y-direction where the am- 
plitude of the eigenfunctions is close to zero. We used a 
grid size of 160 x 120 for the Gaussian integration. This 
grid size is sufficiently large such that the numerical er- 
ror in the Gaussian quadrature is much smaller than the 
error due to basis truncations. We note however, that 
for a larger magnetic field, a bigger grid might be re- 
quired due to more rapid fiuctuations in the phase of the 
eigenfunctions V'i(r)- To make sure that the y cutoff is 
reasonable and the grid is sufficiently dense we checked 
the normalization of the eigenfunctions. 



The Coulomb potential ( 10 1 is integrable in the origin, 
i. e. for r = r', in two dimensions, for 77 = 0. Therefore 
the integral ([9| is mathematically convergent and the 
regularization parameter rj is theoretically not needed. 
However, due to the discretization of the two-dimensional 
space, working in practice with 77 = can nevertheless 
cause problems in the numerical integration. A quick 
way around this problem is replacing Ijrij) with X,>(r) 
where 



X,>(r)^ /{V;(r')-V;(r)} 
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r — r'\ + 1] 



(11) 



It's easy to show that the transformation Tjr{r) — > ^jr{r) 
leaves He unchanged and conveniently rids of us of the 
convergence problems we had with Xjr{r). The validity 
of this transformation does not depend on geometry or 
dimension |16| p. 63-65]. We note that even though the 
limit ?7 — >■ 0+ is well defined in pTj ) , we still have to keep 
?7 > for numerical reasons. However, we can have rj 
much smaller than if we used (lol directly. 

Now that we have the form of the many-electron 
Hamiltonian we need to find a suitable basis for the 
many-electron Fock space. The natural choice is the oc- 
cupation number basis {|/i)} where 
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which means that n^ 



particles are in state j^'i 



state |"!/'2) etc. We use Latin indices for the single-electron 



states and Greek ones for the many electron states, 
fermions we have nf = or nf = 1. For example. 



For 
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When doing calculations, the Fock space needs to be 
truncated by putting cx) — >■ A^ses in (12 1, where TVges 
is a finite positive integer. This means we are using a 
finite number of single-electron states to construct the 
Fock space. This is the first truncation we perform on 
Fock space. The corresponding number of many-electron 



states N,r,r, is 






where Np is the number of electrons. 



This rapid growth of N^^^s obviously limits us to calcu- 
lations for a few electrons only. 

To use this Fock basis we need some way to uniquely 
number the states. We need some mapping F : \fi) — >■ /i 
where /i S Z+ and it's inverse F"^ : /i — > |/i). There 
are many ways to construct F. The exact details will 
depend on factors such as whether or not all the states 
l/i) contain the same number of electrons. For a closed 
system the electron number is constant [11^ , however an 
open system would have a varying number of electrons 
|12] . For this reason we will not go into details of the 
form of F, but assume that we have such a mapping. 

We can now calculate the matrix representation of He 
in the {\fi)} basis using 
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where {^\dldJ,dsdr\i') is calculated using jTS] 
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The phase factor (—1)''''' ensures that d] and di satisfy 
the fermionic anti-commutation relations. Next we diag- 
onalize He and find its eigenstates |^) and energies i5^. 
In the diagonalization process we obtain a unitary trans- 
formation V which satisfies 



JVir 



Im) = E ^M^i 



(18) 



We distinguish between the many-body noninteracting 
and the many-body interacting states by using an angular 
bracket for the kets of the first type, |/Lt), and a rounded 
bracket for the kets of the second type, |/x), respectively. 



This unitary transformation will be used extensively 
because it is much more efficient to perform calculations 
in the {\n)} basis and perform a unitary transformation 
to {|/i)}, rather than explicitly calculating and storing 
the many-electron eigenfunctions. This means that ev- 
ery time we need |/i) for calculations, we need to perform 
a unitary transformation using a matrix that has the di- 
mension A'mes X -^mcs ■ This cau be a problem since iVmcs 
is a rapidly increasing function of Ne and A^sos- For our 
calculations we use A^sos — 50 for two electrons, result- 
ing in Amines = {^2) ^ 1225. For three electrons we use 
TVses ^ 30, resulting in iV^es = (^3°) = 4060. The case for 
a single-electron is trivial since Nses — Mnes- For these 
values of A^ses and electron numbers we get a truncation 
error that is smaller than the error due to the trunca- 
tion of the electron-photon Fock space which is covered 
in the next section. For this reason we will not go into 
discussion of convergence for the purely electronic Fock 
space. 

Before we go on and include interaction with a quan- 
tized EM field we note that if two Fock states |^) and 
\v) do not have the same number of electrons, then 
{lJ,\dldjdsdr\i') = {fi\dldj\iy) = for all i,j,r,s. In other 
words the Coulomb interaction conserves the number of 
electrons. This means that there exists a basis where 
He is block diagonal, where each block consists of states 
with the same number of electrons. Therefore, there exist 
unitary transformations Vn^ for each number of electrons 
which has the same dimension as the block of 'He corre- 
sponding to A''e electrons. We can therefore use many 
small unitary transformations for each electron number 
instead of a big one which works for all number of elec- 
trons. This can be a big boost in computation speed for 
large matrices. 



where tt = p + qA is the mechanical momentum. The 



term in ( 20 ) is the paramagnetic interaction term and 



(21 1 is the diamagnetic term. To go further we need 
to decide upon the form of Aem- We assume that the 
single-mode photon wavelength is much larger than char- 
acteristic length scales of the system. We can then ap- 
proximate the vector potential amplitude to be constant 
over the electronic system. Although related, this is not 
exactly the dipole approximation since we will not omit 
the diamagnetic electron-photon interaction term. We 
can then write the vector potential as 

Aem =i eAEM(a + a'^) = e-— ^ — (a-|-a^), (22) 

where e is a unit vector in the direction of the field po- 
larization and £(. = qA-^-^Vt^a-^ is the electron-photon 
coupling strength. 

The strength of the photon-electron coupling is char- 
acterized by Aem, the magnitude of which depends on 
the experimental setup. For a 3D Fabry Perot cavity we 
would have Aem — •\//i/(2wpV^eo) where V is the cavity 
volume. Another potential setup is a ID transmission 
line resonator [5 where it would be more appropriate to 
write Aem in terms of the electric field vacuum fluctua- 
tion El"^^ = v/(0|E-E|0> where E = -dA^M/dt and |0) 
is the lowest eigenstates of "Hem- We would then have 
Aem = El^^ /ujp. 

Using the approximation in Eq. ( 22 I , the expressions 



/(1'2) • 



for 'Hi,-,^ in Eqs. (20 1-(21 1 can be greatly simplified since 



we can pull Aem in front of the integrals and the com- 
mutator [Aem,'''] is zero. For the paramagnetic term, 
we get 



n\2^£c{a + a^)Y.g^Ad, 



(23) 



IV. INCLUSION OF A QUANTIZED EM FIELD 

Now suppose the system described in section |TT] is sub- 
ject to a single-mode quantized electromagnetic field with 
vector potential Aem- We can write the Hamiltonian as 



ri — ric ^" ri 



EM ■ 



Hi: 



(19) 



where T-Lp is the purely electronic Hamiltonian including 
the Coulomb interaction, Hem is the free field photon 
term and Hint contains the electron-photon interaction. 
Ignoring the zero point energy, the free field term can be 
written as Hem = hujpa^a where huip is the single photon 
energy and a (a^) is a bosonic annihilation (creation) 
operator. The electron-photon interaction term can be 
split into two terms Hint = Hij^j 



hU where 



m^l ^ E (^« I ^ (^ • Aem + Aem • tt) |^,)4d, (20) 



^int^E(V'^l|;:IAEMriV',)dk- 
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where gij is the dimensionless coupling between the elec- 
trons and the cavity mode defined by 
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dr[^*(r){7r^,(r)} + {7r^*(r)}^,(r)] 



(24) 



The dimensionless coupling gij is closely related to the 
dipole transition moment dij = ~q{i\r\j) according to 
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A very accurate way to compute gij is to calculate the 
integral in ( |24[ ) analytically in the original one electron 
basis {\(f>n) ^Ifm)} and perform a unitary transforma- 
tion into the {|i/'i)} basis. Another simpler method is to 
store the x and y derivatives of ipi (r) on a grid and cal- 
culate (24) using Gaussian quadrature. This method is 



less accurate but is easier to implement. 
As for the diamagnetic term, we get 



H 



(2) 
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U^ , (26) 



where N'^ is the number operator in the electron Fock 

space. An interesting aspect of "Hj^t i^ thai it contains 
no dependence on the photon polarization or geometry 
of the system. 

A natural choice of basis for doing calculations is 
{|/i) (X) |A/)} = {|a)} where \M) are eigenstates of the 
photon number operator a) a, with M the number of pho- 
tons. We will obviously need another bijection to label 
the states |/i) ® \M) with a single index a. The depen- 
dence of /i and M on a is suppressed for easier read- 
ing. For the {|a)} basis we use the lowest A^mesT <C 
Ames Coulomb interacting eigenstates and photon states 
containing up to A'em photons, resulting in a total of 
A^mesT X (A^EM + 1) statcs in the {\a)} basis. This is the 
second time we truncate a many-body Fock space. Ap- 
propriate values of A^mesT and Aem are investigated in 
section |Vl 

Calculating matrix elements of "He and "Hem is 
straightforward in the {\a)} basis. We get 

(/x; M\'H,\y; N) = E^S^Jmn (27) 

(/i; M|-Hem|j^; N) = MfuvpS^^SMN , (28) 

where the shorthand |^;M) = |/i) (g) \M) has been used. 
For the paramagnetic interaction term we get 

{ti;M\n[llW;N} = £,Y.9^AMdJW){M\a + a^\N} 

ij 

= £cQtiy [yM + 15n,m+i + VA + 15m,n+i] , 

(29) 

where we define 

g^, ^ ^5.,(A*l4rf.» =Y.g^M'^^4d,V\v) , (30) 

ij ij 

which is the many-electron generalization of gij . Its con- 
nection to the electron-photon coupling energy in the 
Jaynes-Cummings model is explained in Ref. [TT]. We 
will refer to Q^i, as the dimensionless geometric coupling 
(DGC) between the electronic states |/i) and \v). As for 
the diamagnetic term we get 



{^,-M\H^\v-N)^^N,5,, 



[N+-)5mn+ 



^ v/(Af + l)(M + 2)<5w.M+2 
i .J{N + 1){N + 2)5m,n+2 



(31) 



where A^ is the number of electrons in the state |/i). The 
matrix elements of the total Hamiltonian (/i; M\'H\v, A) 
are obtained by adding ( 27 1 , ( 28 1 , ( 29 1 and ( 31 ) together. 



The final step is diagonalizing "H and obtaining the 
allowed energies E^ and the corresponding eigenstates 
I a) which are related to \a) by the unitary transformation 



that is obtained in the diagonalization process. Again we 
use the right angular bracket for the basis states and the 
rounded bracket for the interacting states, this time the 
interaction being between electrons and photons. 

Expectation values of an observable A can then be 
calculated using 

(^> = Tr(p^) = ^(a|^|/3)p^„ = ^(a|^|/3)p^„ (33) 

a/3 aP 

where p (p) is the density matrix of the system in the \a) 
{\a)) basis. The main advantage of working in the \a) 
basis is that (Qf|yl|/3) is easy to calculate. However it is 
very hard to truncate p effectively. Working in the \a) 
basis is the exact opposite, (a|^|^) — {a\yV^ AyV\$) is 
expensive to calculate but it's easy to truncate p because 
if the system is in an energetically low state, all of its 
biggest elements are concentrated in its top left corner 
(low a and /3). Note that, although (a|^|^) is relatively 
expensive to calculate, it can be computed beforehand 
and saved. 

Example of an interesting observable is the photon 
number operator A/'^'^ — a^a. Its expectation value can 
be calculated using 



{MP'')^J2{p;M\a^ay,N)pp^ 
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(34) 



Another interesting observable is the charge density 

Q(r) = -g^V*(r)V',(r)dk-, (35) 



/3 



(32) 



the expectation value of which can be calculated using 

(2) ('^) = EE^*WV',(r)(MMlrf.»'5AfArPa0 , (36) 
ap ij 

where it is important to calculate {p,\d\dj\v) = 
{pIV^ dldjV\i^) beforehand to avoid unnecessary repeti- 
tions. 



V. RESULTS 

For the results presented in this section we use B = 
0.1 T, hno = 1.0 meV, L^ = 300 nm, m = 0.067TOe 
and e = 12.4eo (GaAs parameters). We choose cOp such 
that the system is on resonance between some chosen 
electron states \k) and |A) with detuning 6, that is hutp — 
\Ex -E,,\+S where 6 = 0.01(i?A - E^). We refer to \k) 
and |A) as the active states. We use A > k, making S 
positive so that we are slightly over resonance. Choosing 
A < K would give a negative 6, resulting in a system 
that is slightly under resonance. To distinguish between 
electron states with different number of electrons, we use 



the notation \^)n^ to denote the /i-th electronic state 
containing Ng electrons. For example. |4)2 is the fourth 
lowest two electron state. 

Figure [T] shows the energy spectra of H as a function 
of the coupling strength £c for both x and y polarization. 
The importance of the diamagnetic interaction term is 
also illustrated in Figure IT] by plotting the same en- 
ergy spectrum, but omitting the diamagnetic term. For 
small coupling, ignoring the diamagnetic term is a valid 
approximation. However, for higher coupling strength, 
the model without the diamagnetic term starts exhibit- 
ing red shift with respect to the exact result. This red 
shift becomes visible at around I^^aI £-c/^p ~ 0.1, where 
the values of A, k, \Qk,\\ and fujp are given in the figure 
text. For even higher coupling strength, the results with- 
out the diamagnetic term start exhibiting an unphysical 
downwards dive in energy. In this regime, the results 
are highly divergent with respect to iVmosT- However, 
keeping A^mosT constant, the results are convergent with 
respect to TVem- 

To get an estimate of the numerical truncation errors 
we look at the relative variation of the energy of state |a) 
defined as 



a) 



K 



(a) _ 



r(") _ £;(") 



E. 



(a) 



(37) 
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where E^"' is the energy of state \a) and i refers to a spe- 
cific parameter related to the size of the truncated Fock 
space. For example i can be N^^s, -^mcsT or N-em- Typ- 
ically, j is the maximum value of that parameter which 
can be used to obtain the numerical output in a reason- 
able computing time. We vary i and check the converge 
of the results. When changing the parameters i and j, 
all other accuracy parameters are kept constant. We also 
define the maximum error of the N lowest states as 



pmax 



/96[1,JV] 






(38) 



The value we choose for TV depends on what we intend 
to use the states for, once we have obtained them. For 
calculating electron transport using the generalized mas- 
ter equation 64 states are typically used, so that is the 
value we will use for N |12| . Our criteria for convergent 
results is that the error is not visible on a graph such as 
in Figure [l] This condition translates into a maximum 
relative error of ~ 10""^. 

Figure |2] shows the relative error in the energy spec- 
trum for one electron due to the finite value of A^mosT, 
that is the error due to the truncation of the electron 
part of the joint electron-photon Fock space basis. From 
the figure we see that for A^mesT = 200, results are con- 
vergent up to \Gi2\£c/f>^p — 0.6 for a;-polarization and 
IGi^l £c/fuxjp — 0.7 for y-polarization. From the figure we 
also see that the error rises very rapidly for small £c but 
as £c becomes a considerable fraction of hujp, the error 
increases much slower. 

Figure [3] shows the relative error due to the finite 
value of A'mesT for two electrons and both polarizations. 
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Figure 1. Energy spectra for the lowest ~ 60 states with 
one electron and a;-polarization (a) and two electrons and y- 
polarization (b). The diamagnetic A'^ term in the e-EM inter- 
action Hamiltonian is both included (blue) and omitted (red). 
In (a), the system is on resonance between the one electron 
states |l)i and |2)i with a DGC strength of |5i2| = 0.290 and 
hujp — 0.185 meV. In (b), the system is on resonance between 
the two electron states |1)2 and |5)2 with a DGC strength of 
tjisl — 0.987 and hcOp = 1.025 meV. As can be seen from 
the figure, omitting the A^ term does give accurate results 
for small £c, while for large £c the energy spectrum takes a 
steep dive downwards. This dive also takes place in the two 
electron case, however it can't be seen in the chosen range of 
£c. There is no physical significance in these dives since the 
results are highly divergent in those areas. 



For A^mcsT = 200, the results are convergent up to 
\Gi2\£c/fi^p — 0.3 for x-polarization and \Gi5\£c/fi^p — 
0.3 for y-polarization. The same convergence calcula- 
tions for 3 electrons (not shown here) gives convergent 
results for \Gi3\£c/fiijJp — 0.25 for cc-polarization and 
1^15 1 £c/f^p — 0.25 for y-polarization. 

Figure |4] shows the relative error due to the finite value 
of A^EMi that is the truncation of the photon part of the 
joint electron-photon Fock space basis. From the figure 
we see that a modest value of TVem = 20 is enough for 
the error to be 3 — 12 orders of magnitude smaller than 
the A'liicsT truncation error shown in figures [2] and [3] The 
results in Figure [4] are for one electron but the two and 
three electron cases (not shown here) exhibit the same be- 
havior. The reason for this faster convergence w.r.t A^em 
is most likely that the electronic energy spectrum is much 
more dense, with a high amount of energy crossings/anti- 
crossings, which requires a larger basis. 
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Figure 2. Convergence calculations with respect to NmcsT for 
a;-polarization (a) and {/-polarization (b). In (a), the system 
is on resonance between the one electron states |l)i and |2)i 
giving hiUp = 0.185 meV and \Qi2\ ~ 0.290. The results are 
convergent up to £c — 0.37 or |t/i2| £c/fiijJp — 0.6. In (b) the 
system is on resonance between the one electron states |l)i 
and |5)i giving hujp = 1.03 meV and |5i5| =0.701. The results 
are convergent up to £c — 1.05 or \Gi5\ Ec/hup ~ 0.7. For this 
run we have a — 100, b — 150, c = 200 and d — 250 (see 
equations |37] and [38] for definition) . The maximum number 
of photons is kept constant at TVem = 20. 
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Figure 3. Convergence calculations with respect to A'mcsT for 
x-polarization (a) and j/-polarization (b). In (a), the system 
is on resonance between the two electron states |1)2 and |2)2 
giving hiOp = 0.516 meV and \Q\2\ ~ 0.648. The results are 
convergent up to £c — 0.25 or \Qi2\£c/fvjjp ~ 0.3. In (b) 
the system is on resonance between the two electron states 
11)2 and |5)2 giving hcup = 1.03 meV and l^isl = 0.987. The 
results are convergent up to £c — 0.25 or jSisI £c/huip — 0.24. 
For this run we have a — 100, b = 150, c = 200 and d = 



250 (see equations 37 and |38| for definition). Other accuracy 
parameters are Nscs = 50 and TVem = 20. 
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Figure 4. Convergence calculations with respect to A'^em for 
x-polarization (a) and y-polarization (b). Values of hujp and 
\G\k\ are the same as in Figure [2] for both polarizations. We 
can see that for Nem = 20 (green) , the results are acceptable 
for the whole range of £c considered. For this run we have 
a = 10, 6 = 15, c = 20 and d = 25 (see equations 37 and 38 
for definition) . The electron state number is kept constant at 

NmcsT = 200. 
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Figure 5. Convergence calculations for two electrons with 
respect to A'^mesT for x-polarization (a) and j/-polarization (b) . 
For both polarizations, the system is off resonance with hwp = 
0.4 meV. In both cases, the results are convergent up to fc — 
0.26 or £c/hujp ~ 0.65. For this run we have a = 100, b = 150, 



c = 200 and d = 250 (see equations 37 and 38 for definition) 
Other accuracy parameters are A'ses = 50 and A'^em = 20. 



Although in this paper we have put the photon fre- 
quency on resonance between two electronic states, we 
are in no way forced to do so (see Ref. [H]). This moti- 
vates us to investigate convergence for a system that is off 
resonance. Figure [s] shows convergence calculations for a 
system that is off resonance and contains two electrons. 
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From the figure we see that the results are convergent up 
to Ec/fiiiJp — 0.65 for both x and y polarizations. The 
reason we use the ratio Ed^-p rather than \Q^\\ Ec/fu^jp 
is that when the system is ofi^ resonance, the concept of 
active states |A) and \k) has no meaning. 



VI. CONCLUDING REMARKS 

We have described a rigorous method to compute the 
many-body states of a multi level Coulomb interacting 
electronic system which also interacts with a single-mode 
quantized EM field. The model is exact in the sense 
that the only approximations are the finite size of the 
single- and many-body bases and the finite size of grids 
on which single-electron eigenfunctions are stored. The 
convergence with respect to these parameters is carefully 
controlled. 

Due to the exact numerical nature of the model, cal- 
culations for arbitrarily strong photon-matter interac- 
tion can in principle be performed with a big enough 
basis. Numerical results show that the main bottle- 
neck is the large number of electron states needed in 
the joint photon-electron many-body basis. Convergence 



with respect to the number of photon states is much 
faster where ^ 20 states are sufficient to guarantee nu- 
merical error that is 3 — 12 orders of magnitude smaller 
than the error caused by the electronic basis truncation 
with ~ 200 states. We have found that including the 
diamagnetic photon-electron interaction term drastically 
improves convergence when the electron-photon coupling 
strength is considerable in size to the single photon en- 
ergy (ultrastrong coupling regime). Without the diamag- 
netic term, the model shows unphysical behavior in the 
ultrastrong coupling regime due to divergent results. 
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